Confirming and extending the hypothesis of universality in sandpiles 
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Stochastic sandpiles self-organize to a critical point with scaling behavior different from directed 
percolation (DP) and characterized by the presence of an additional conservation law. This is usually 
called C-DP or Manna universality class. There remains, however, an exception to this universality 
principle: a sandpile automaton introduced by Maslov and Zhang, which was claimed to be in the 
directed percolation class despite of the existence of a conservation law. In this paper we show, 
by means of careful numerical simulations as well as by constructing and analyzing a field theory, 
that (contrarily to what previously thought) this sandpile is also in the C-DP or Manna class. This 
confirms the hypothesis of universality for stochastic sandpiles, and gives rise to a fully coherent 
picture of self-organized criticality in systems with a conservation law. In passing, we obtain a 
number of results for the C-DP class and introduce a new strategy to easily discriminate between 
DP and C-DP scaling. 

PACS numbers: 05.50.+q,02.50.-r,64.60.Ht,05.70.Ln 



I. INTRODUCTION 

Aimed at shedding some light on the origin of scale 
invariance in many contexts in Nature, different mecha- 
nisms for self-organized criticality (SOC) were proposed 
during the last two decades following the seminal work 
by Per Bak and coauthors P, SHI- Sandpiles, ricepiles, 
and earthquake toy-models become paradigmatic exam- 
ples capturing the essence of self-organization to scale- 
invariant (critical) behavior without apparently requiring 
the fine tuning of parameters (see 0, [a, H] ) • 

Sandpiles played a central role in the development of 
this field [H, H, El • They are metaphors of real systems 
(as earthquakes, snow avalanches, stick-slip phenomena, 
etc.) in which some type of stress or energy is accu- 
mulated at some slow timescale and relaxed in a much 
faster way. In sandpiles, grains are slowly added and 
eventually relaxed whenever a local instability thresh- 
old is overcome; then they are transmitted to neigh- 
boring sites which, on their turn, may become unstable 
and relax, generating avalanches of activity. Considering 
open boundaries to allow for an energy balance, a steady 
state with power-law distributed avalanches is eventually 
reached. 

In order to rationalize sandpiles in particular, and SOC 
in general, and to understand their critical properties, it 
was proposed @, [1, Q to look at them as systems with 
many absorbing states [l(| [HI ■ The underlying idea is 
that, in the absence of external driving, sandpile mod- 
els get eventually trapped into stable configurations from 
which they cannot escape, i.e. absorbing states [To] . The 
concept of fixed energy sandpiles (FES) was introduced 
to make this connection more explicit. FES sandpiles 
share the microscopic rules with their standard (slowly 
driven and dissipative) counterparts, but with no driving 
nor dissipation. In this way, the total amount of sand 
or energy becomes a conserved quantity acting as a con- 



trol parameter. Calling E c the average density of energy 
of a standard sandpile in its stationary (self-organized) 
critical state, it has been shown that the corresponding 
fixed-energy sandpile exhibits a transition from an active 
phase to an absorbing one at, precisely, E c . It can be 
argued that slow driving and dissipation acting together 
at infinitely separated time-scales constitute a mechanism 
able to pin a generic system with absorbing states and a 
conservation law to its critical point [3, 0, @, H, fl2| . 

Usingthe relation with transitions into absorbing 
states [13[ , a Langevin equation describing FES stochas- 
tic sandpiles was proposed [1, Q in the spirit of Ho- 
henberg and Halperin [Til ]. The Langevin equation 
is similar to the well known directed-percolation (DP) 
Langevin equation describing generic systems with ab- 
sorbing states [13, EH but, additionally, it is coupled 
linearly to a conserved non-diffusive energy field, repre- 
senting the conservation of energy grains in the sandpile 
dynamics [!,[§]: 

dtp(x,t) = ap — bp 2 + ujpE(x,t) + V 2 p + a^/prj(x,t) 
d t E(x,t) = DV 2 p, (1) 

where D, a, b and u> are constants, p(x, t) and E(x, t) 
are the activity and the energy field, respectively, while 
r\ is a zero-mean Gaussian white noise. The universal- 
ity class described by this Langevin equation, which in- 
cludes sandpiles as well as all other systems with many 
absorbing states and a non-diffusive conserved field, is 
usually called Manna or C-DP class [1 i, [II E3- As 
a side note, let us remark that the original Bak-Tang- 
Wiesenfeld model, being deterministic, has many other 
conservation laws (toppling invariants) and is therefore 
not described by the present stochastic theory [H, [TU . 

Incidentally, even if it has been clearly established that 
DP and C-DP constitute two different universality classes 
[l9l . I20I [2l| , most of their universal features (critical ex- 
ponents, moment ratios, scaling functions,...) are very 
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similar, making it difficult to discriminate numerically 
between both classes in any spatial dimension. 

Despite of the fact that the critical behavior of stochas- 
tic sandpiles is accepted to be universal and described by 
Eq.([T]), there are a few sandpile models which, despite of 
being conservative, have been argued to exhibit a differ- 
ent type of critical behavior [12, namely DP, violating 
the conjecture of universality: 

• the Mohanty-Dhar (MD) sandpile 22], in which 
grains have some probability to remain stable even 
if they are above the instability threshold. It was 
argued, based on numerical simulations and on ex- 
act results for an anisotropic version of it, that the 
(isotropic) MD sandpile should exhibit DP-like be- 
havior [22| 

• the Maslov and Zhang (MZ) sandpile (23j . in which 
not only the grains of unstable sites are redis- 
tributed among nearest neighbors. Instead, all the 
grains in the local neighborhood of unstable sites 
are randomly redistributed or reshuffled. On the 
basis of Monte Carlo simulations this sandpile was 
argued to exhibit DP behavior. 

Aimed at clarifying this puzzling situation, in a recent 
work we provided strong numerical and analytical evi- 
dence that, contrarily to what previously thought, the 
(isotropic) MD sandpile is actually in the C-DP class 
[24],[25[ (see [26[ for a discrepant viewpoint). To further 
substantiate our claim, in [2|| we proposed a strategy 
to easily discriminate between DP and C-DP consisting 
in introducing a wall (either absorbing or reflecting); sys- 
tems in the DP and in the C-DP classes behave in (quali- 
tative and quantitatively) very different ways in the pres- 
ence of walls, providing an easy criterion to discriminate 
between both classes. 

Our goal in the present paper is to scrutinize the last 
remaining piece in the puzzle, i.e. the MZ model, by em- 
ploying extensive numerical simulations as well as field 
theoretical considerations. For the numerics, we take ad- 
vantage of the previously introduced discrimination crite- 
rion and, also, present another method to easily discrim- 
inate between DP and C-DP, based in the introduction 
of anisotropy in one spatial direction. 

The analyses presented on what follows, show in a 
clean-cut way that the MZ model is actually in the C- 
DP class. This leaves no dangling end in the sandpile 
universality picture and, as a byproduct, confirms the 
general validity of the absorbing-state approach to ra- 
tionalize SOC. New results for the C-DP class are also 
obtained. 



II. THE MASLOV-ZHANG SANDPILE 

The Maslov-Zhang (MZ) sandpile [23[ is defined as fol- 
lows: 



1. Driving: an input energy, 5E < 1, is added to the 
central (or to a randomly chosen) site, i, of a d- 
dimensional lattice, and the site is declared active. 

2. Relaxation: energy is locally redistributed ( "reshuf- 
fled") between the active site and its nearest- 
neighbors, according to: 

r 2d+1 

E i = ^2d+l H ( 2 ) 

where are uniformly distributed (r, € [0, 1]) ran- 
dom variables, and the sums are performed over the 
site i and its 2d nearest neighbors. 

3. Activation: each of the sites involved in the reshuf- 
fling is declared active with a probability given 
by its own energy, triggering the generation of 
avalanches of activity. 

4. Avalanches: new active sites are added to a list and 
relaxed in a sequential way. 

5. Dissipation: energy arriving at the open borders is 
removed from the system. 

6. Avalanches proceed until all activity ceases, and 
then a new external input is added. Eventually a 
critical stationary state is reached. 

Monte Carlo simulations by Maslov and Zhang re- 
vealed exponents compatible with those of directed per- 
colation in two dimensions and above, while in one di- 
mension some anomalies in the scaling were reported 
[23l ]. Later, simulations of the FES counterpart of the 
MZ sandpile led to very similar results Q. 

Before entering a more careful numerical analysis of the 
MZ sandpile, we take a detour to construct explicitely a 
Langevin equation for such a model. This will allow us 
to better understand what the main relevant ingredients 
of the MZ model are and in what sense it differs from 
other C-DP models and from Eq.|T|). 



III. A LANGEVIN EQUATION FOR THE 
MASLOV-ZHANG SANDPILE 

The main difference between the Maslov-Zhang cellu- 
lar automaton and other sandpiles in the C-DP univer- 
sality class is that, while in the C-DP class as the Manna 
model, the only energy redistributed by the dynamics (by 
topplings) is that accumulated in active sites, in the MZ 
dynamics both the energy of active sites as well as that 
of its nearest neighbors is redistributed. This leads to a 
more severe local redistribution of energy, that we quan- 
tify in what follows in terms of a new set of Langevin 
equations, first in a phenomenological way and then by 
deriving it from the microscopic rules. 
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A. Phenomenological Langevin equation 

Diffusion of energy occurs in the C-DP class by means 
of activity relaxation. At a mesoscopic level, this implies 
that the rate of change of the energy density E(x, t) at a 
given position is proportional to minus the divergence of 
a current, d t E(x,t) — —W.j(x,t), where j(x,t) is given 
by the gradient of the activity field, 



j(x,t) = -DVp(x,t) 



(3) 



leading to d t E(x, t) = DV 2 p(x, t) in Eq.©. 

Instead, in the MZ model, changes of energy are con- 
trolled by the reshuffling rule, i.e. energy does not need 
to be at an active site (but in its neighborhood) to be 
redistributed. Hence, at a mesoscopic level the current 
is given by gradients of the energy in the presence of 
non-vanishing activity. More specifically, dtE(x, t) = 
— V.j(x, t) where now 



j{x,t)=-DVE(x,t), 



(4) 



and the diffusion D is not a constant but a functional pro- 
portional to p(x,t), D(p(x,i)) = Dp(x,t) capturing the 
requirement that local reshuffling of energy only occurs 
in the presence of activity. This enforces the absorbing 
state condition that dynamics ceases if p = 0. This finally 
leads to 



d t E(x,t) = D V.[p(x,t)VE(x,t)]. 



(5) 



On the other hand, the equation for the activity is not 
expected to change in any relevant way, so one obtains 

dtp{x, t) = ap — bp 2 + V 2 p + u>pE(x, t) + a^/prj(x, i) 
d t E{x,t) = DV.\p{x,t)VE{x,t)] (6) 

representing the MZ dynamics at a mesoscopic level. 
This is to be compared with Eq. (fTJ) . 



B. Microscopic derivation of the Langevin equation 

To gain more confidence on the phenomenological set 
of equations, Eq. ([6|) let us now present an explicit mi- 
croscopic derivation. We consider a parallel version of 
the MZ model in which all active sites are relaxed at ev- 
ery time step. To do so, we assume that each site uses a 
fraction 1 / (2d+l) of its total energy for eventual redistri- 
butions with each of the sites in its local neighborhood. 
At mean-field level, the energy evolves according to 



2d+l 
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i,t + l 



2d 
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where j (respectively k) runs over the nearest neighbors 
of i (j) and 9 is the Heaviside step function (0(z) = 



if z < 0) . For any inactive site in the local neighborhood 
(i.e. pj = 0), the central site does not redistribute the 
corresponding fraction of its energy (first term in Eq.©). 
For any active site in the local neighborhood (i.e. pj > 
0) the central site receives on average a corresponding 
fraction of energy (second term). Eq.© is a mean-field 
equation to which fluctuations should be added to have 
a more detailed description. Such fluctuations appear 
as a (conserved) noise for the resulting energy equation 
and can be easily argued to constitute a higher order, 
irrelevant, correction. 

Regularizing the Heaviside step function in Eq.© by 
means of a hyperbolic-tangent and expanding it in power 
series up to first order in p [24j], we obtain 
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3=1 



Introducing the d-dimensional discrete Laplacian ~S/ 2 Ei 
yjj—i {Ej — Ei) the first term can be rewritten as: 



2d 
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2d+l 

Similarly, the second term can be expressed as 
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(2d+l) 



ftV 2 ^+ (2rf + 1)2 v 2 ( Pl v 2 ^). 

Putting these two contributions together and reorganiz- 
ing the discrete derivatives, one obtains: 
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which in the continuous time limit becomes 



(11) 



d t E t 



2d+l 



V. {p l ^E l , 



{ 2d + —^ 2 ^ 2 ^)- 
(12) 

The second term in this last equation can be argued to be 
naively irrelevant in the renormalization group sense (as 
it includes higher order derivatives) and hence dropped 
out. Identifying D with 2/(2d+ 1) we recover the phe- 
nomenological Eq.([5]), as the leading contribution. 

As said above, this is to be compared with the stan- 
dard Langevin Eq.([T]) for the C-DP class, in which the 
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flux of energy is given by the gradient of the activity 
itself. The following questions pop up naturally: Does 
Eq. (|6|) lead to a critical behavior different from that of 
Eq. (TT])? Could this be the reason why the MZ sandpilc 
was claimed to differ from the C-DP class? Is this new 
form of conserved energy dynamics irrelevant at the DP 
fixed point, supporting the MZ sandpile to be DP like? 

To properly answer these questions one should resort 
to a full renormalization group calculation. Given that 
even for the C-DP class this has proven to be a, still 
not satisfactorily accomplished, difficult task [13] , we will 
leave aside such a strategy here. Instead, in section IV 
we will give an answer to these question by means of 
computational studies of the microscopic MZ model as 
well as of its equivalent Langevin Eq.©. 



C. Naive power counting 

A power counting analysis is not helpful in elucidating 
the relevancy of the energy-diffusion term, Eq. ([5]) . Actu- 
ally, as there is a lineal dependence on the energy field in 
both the r.h.s. and the l.h.s. of Eq.fTj, there is no way 
to extract the energy field naive dimension nor, hence, to 
make any statement about the relevancy of the coupling 
term wp(x,t)E(x,t) at the DP fixed point. 

On the other hand, it is easy to cast Eq.© into a 
generating functional following standard procedures, to 
set the basis of a perturbative expansion. However, one 
soon realizes that technical difficulties similar to those 
encountered for the analysis of Eq. © (including the pres- 
ence of generically singular propagators [1,113) snow U P) 
hindering the perturbative calculation. It is our believe 
that some type of non-perturbative technique, or non- 
conventional perturbative expansion, is required to elu- 
cidate the renormalization group fixed point of this type 
of problems. 

The analytical understanding, at a field theory level of 
C-DP as well as the MZ-Langevin equation remain open 
challenging tasks. 



IV. MONTE CARLO SIMULATIONS OF THE 
MASLOV-ZHANG SANDPILE 

We have performed extensive Monte Carlo simulations 
of the MZ sandpile (and variations of it) and scrutinized 
its asymptotic (long time and large system size) prop- 
erties. We report on two different types of numerical 
experiments. 

• "SOC" or avalanche experiments: 

By iterating slow addition of grains in the sandpile 
with open boundaries, the system self-organizes 
to a state with average energy E c . Then, the 
avalanche size distribution, P(s) and the avalanche 
time distribution -Pt(t), can be estimated and their 
corresponding exponents, r and Tt [28| , measured. 



• Absorbing state experiments: 

At the stationary state, we perform spreading ex- 
periments from a localized seed, and we measure: 
i) the mean quadratic distance to the initial seed 
R 2 r~j t Zapr in active runs, ii) the average number of 
active sites as a function of time, N(t) ~ t v and iii) 
the surviving probability up to time t, P s ~ t~ s . 
We also study the decay at criticality of a homo- 
geneous initial activity (p(t) ~ t~ e ) in the fixed 
energy case [l(| HH ■ 

Scaling laws relating avalanche to spreading exponents 
where described systematically in [28j : two of them are: 
r = (l + r] + 2S)/(l + r] + S),Tt = 1 + S. We measure the 
exponents independently and use scaling laws as a check 
for consistency. 

We simulate the MZ automaton in one-dimensional lat- 
tices up to size L = 2 15 . The stationary critical energy 
density is EL = 0.4928(2) and, contrarily to what re- 
ported in [23[, we do observe clean scaling at criticality, 
even if it emerges only after significatively long transients 
(results not shown) justifying why smaller scale simula- 
tions can lead to erroneous conclusions. 

The resulting critical exponents are gathered together 
in TableQ] They are closer in all cases to the C-DP values 
than to DP ones. Still, given the similarity between the 
numerical values in both classes, it is not save to extract a 
definitive conclusion from these results. Higher numerical 
precision would be required to produce fully convincing 
evidence. As said above, the numerical values of DP and 
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DP 


0.313(1) 


0.159(1) 


1.108(1) 


1.159(1) 


1.265(1) 


0.159(1) 


C-DP 


0.350(5) 


0.170(5) 


1.11(2) 


1.17(2) 


1.39(1) 


0.125(1) 


MZ 


0.32(5) 


0.20(5) 


1.13(5) 


1.20(5) 


1.40(5) 


0.13(1) 



TABLE I: Critical exponents for DP, C-DP, and the MZ sandpile 
in one dimension. DP and C-DP values taken from 20, 25, 281. 



C-DP exponents are closer and closer as the dimensional- 
ity is increased (they coincide above d c = 4). Therefore, 
performing numerical simulations, as the ones above, to 
discriminate between both classes by using larger and 
larger times and system sizes in d > 2 , is not a clever 
idea. Instead, it is advisable to use/devise more effective 
numerical strategies to discriminate between DP and C- 
DP in a simple, efficient, and numerically inexpensive 
way. For this we use: 

• the method devised in [25| consisting in analyzing 
how the system responds to the presence of a wall, 
or 

• a new strategy, which exploits the fact that systems 
in these two classes react in remarkably different 
ways to the introduction of anisotropy in space. 

Both of these strategies allow to obtain clean-cut re- 
sults, as shown in the forthcoming two subsections. 
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A. Boundary Driven Experiments 

The influence of walls in systems in the DP class has 
been profusely analyzed in the literature (29J. In par- 
ticular, it is well known that, if spreading (and SOC) 
experiments are performed nearby a wall, the surviving 
probability is significatively affected, and avalanche and 
spreading exponents change in a non-trivial way with re- 
spect to their bulk counterparts [291 ] . It is also known 
that in the DP class, both reflecting and absorbing walls 
lead to a common type of universal "surface critical be- 
havior", that we call surface directed percolation (SDP), 
characterized (in one-dimension) by the exponents shown 
in Table II. 

In contrast, the effect of walls in C-DP systems has 
been studied only recently [25j| . Contrarily to the DP 
case, absorbing and reflecting walls induce different types 
of surface critical behavior. As illustrated in Table II, all 
spreading and avalanche exponents take distinct values 
for an absorbing and for a reflecting wall. Furthermore, 
the numerical differences between the exponents for ei- 
ther type of wall with respect to their corresponding SDP 
counterparts are very large, allowing for easy numerical 
discrimination [2a |. Finally, in the C-DP class, the ex- 
ponents in the presence of a reflecting wall coincide with 
their bulk counterparts 25]. These features imply that, 
by introducing a wall in a given system with absorbing 
states, it becomes straightforward to distinguish if it is 
in the DP or in the C-DP class, with moderate compu- 
tational cost. 

Following this strategy, we simulated the one- 
dimensional MZ sandpile, as defined above, in the pres- 
ence of both reflecting and absorbing walls. In both cases 
a wall is introduced at the origin (position i = 0), and the 
sandpile is studied in the positive half lattice. In the re- 
flecting case, the energy that should go after reshuffling, 
to the leftmost site, at i = 0, (whose energy is fixed to 
zero) is instead added to its closest nearest neighbor to 
the right, i = 1. On the other hand, the absorbing condi- 
tion is imposed by fixing the energy of the leftmost site 
to zero after every iteration of the microscopic sandpile 
rules, i.e. by removing from the system at every iteration 
all the energy received by the leftmost site. 

Fig-© an d Fig.© show the results of simulations per- 
formed in lattices of system size L — 2 15 , averaging over 
up to 810 7 runs. The corresponding exponents are sum- 
marized in Table II. All of them coincide within numerical 
accuracy with the expected values for the C-DP class in 
the presence of reflecting or absorbing walls respectively, 
and differ significatively from those of SDP. For example, 
in the presence of a reflecting (respectively, absorbing) 
wall, the measured value of rj is 0.37 (—0.37) in good 
agreement with the C-DP expectation, r\ — 0.35 (—0.33) 
and in blatant disagreement with the corresponding DP 
value, 77 « 0.046 (« 0.045), which is one order of mag- 
nitude smaller (and of opposite sign in the case of an 
absorbing wall). Similar large differences are measured 
for all the exponents (see table II). Note also that, as is 
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FIG. 1: Avalanche exponents for the one-dimensional MZ 
sandpile in the presence of a reflecting wall. Up: spreading 
experiments; N(t) is the total number of active sites and P s (t) 
is the surviving probability; system size L — 2 lj number of 
runs 500 . Down: avalanche size (main plot) and time (inset) 
distributions; system size L = 2 number of runs 10 5 . Green 
(red) lines mark DP (C-DP) scaling (see numerical values in 
Table E). 
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0.046(2) 
0.045(2) 


0.425(2) 
0.426(2) 


1.25(3) 
1.28(3) 


1.43(3) 
1.426(2) 


1.257(2) 
1.276(2) 


0.154(3) 
0.165(2) 
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0.35(1) 
-0.33(2) 


0.16(1) 
0.85(2) 


1.11(2) 
1.56(2) 


1.15(2) 
1.81(2) 


1.41(1) 
1.43(2) 


0.12(1) 
0.13(1) 


MZ re/ 
MZ abs 


0.37(5) 
-0.36(5) 


0.15(5) 
0.84(5) 


1.09(5) 
1.57(5) 


1.14(5) 
1.84(5) 


1.30(5) 
1.25(5) 





TABLE II: One-dimensional critical exponents for DP and C-DP 
without walls |20l . |28| and in the presence of absorbing and reflect- 
ing walls [25ll . Values in rows 1 (DP re ^) and 2 (DP a (, s ) coincide 
within error-bars. Note also, that values in row 3 ( C-DP re ^ ) coin- 
cide with those for C-DP (table [0. Results for the MZ sandpile in 
the presence of reflecting/absorbing walls are reported in the last 
two rows, theta has not been measured in these cases. 



the case in the C-DP class [25| , the exponents in the pres- 
ence of a reflecting wall coincide within error-bars with 
their bulk counterparts. In conclusion, studying the in- 
fluence of walls we conclude that the one-dimensional MZ 
sandpile exhibits C-DP scaling. 



6 



10" 



10 



10 



1 1 






-N(t) 




i -Ps(t) 




DP 




- C-DP 









P(s)\ 



10 



10 



10 




t 10 



t io 4 : 



10 



10 



10 



FIG. 2: As in Fig. [TJbut for an absorbing wall (numerical 
values summarized in Table ITT)) . System size L — 2 15 and 
8i0 7 runs. 



B. Anisotropic Experiments 

It is well known that systems in the DP class are invari- 
ant under Galilean transformations: if particles have a 
tendency to move anisotropically in one preferred spatial 
direction, that does not alter the critical properties 
The presence of any degree of anisotropy in DP-like sys- 
tems is an irrelevant trait, or in other words, anisotropic 
DP (A-DP) is just DP. 

The role of anisotropy in sandpiles has also been pro- 
fusely studied after the pioneering exact solution by Dhar 
and Ramaswamy (30j of the totally anisotropic or "di- 
rected" counterpart of the Bak-Tang-Wiescnfcld sand- 
pile. Anisotropic stochastic sandpiles have also been 



studied using general principles [3l| and through in- 
terfacial representations [12]. The conclusion is that 
all anisotropic sandpiles, as long as they are stochastic 
[HI], belong to the same universality class, that we call 
anisotropic C-DP (A-C-DP) [34]. The critical exponents 
of models in this class, where first measured numerically 
[34| . and then exactly calculated in any dimension (35| 
(see table Hill and table H¥|) . 

The strategy to be used is straightforward: take the 
MZ sandpile model and switch on anisotropy; if the 
isotropic model was in the DP class, anisotropy should be 
an irrelevant ingredient and the anisotropic counterpart 
should also be DP like. If, instead, the isotropic model is 
in the C-DP class, then anisotropy is a relevant ingredi- 
ent and critical exponent change from C-DP to A-C-DP 
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0.33(2) 


0.14(2) 


1.09(3) 


1.14(3) 


2.00(2) 


A-C-DP 





1/2 


4/3 


3/2 


2 


A-MZ 


-0.02(3) 


0.51(3) 


1.35(5) 


1.48(5) 


1.98(3) 



TABLE III: One dimensional critical exponents for D P |20| . |28H . 
C-DP with a preferred direction (analytical results from [331), an d 
the anisotropic MZ model. 
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DP 


0.230(1) 


0.451(1) 


1.268(1) 


1.451(1) 


1.13(2) 
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-0.07(5) 


0.75(5) 


1.49(5) 


1.70(5) 


1.93(5) 



TABLE IV: Two-dimensional critical exponents for anisotropic 
models: A-DP (i.e. DP), A-C-DP (analytical results), and 
anisotropic-MZ. 



values. 

The simplest way to define an anisotropic MZ (A-MZ) 
model is by fixing one of the rj in Eq.((2|), say the one to 
the right, to its maximum possible value, rj = X, and let- 
ting the others rj to take randomly distributed values in 
[0, 1]. This generates an overall energy flow towards the 
preferred direction (to the right, in this case). Anisotropy 
can be introduced in other ways, including full-anisotropy 
or directness, but this does not affect our conclusions in 
any significant way. 

Fig. [3] and table Hill show our main results for the one- 
dimensional MZ model with anisotropy. Both avalanche 
and spreading exponents are very different from their 
isotropic counterparts. They also differ notoriously from 
DP values, but coincide within error-bars with the ex- 
pected values for the A-C-DP class. The same conclusion 
holds in two dimensions (see Table HV)) . In this way, as 
the anisotropic MZ model belongs to A-C-DP class the 
original, isotropic, MZ sandpile model can be safely con- 
cluded to be in the C-DP universality class, confirming 
the result above. 



V. NUMERICAL INTEGRATION OF THE MZ 
LANGEVIN EQUATION 

In this section we verify that Langevin Eq.© is a 
sound description of the MZ model and that, despite of 
its different form, it behaves asymptotically as Eq.([T]). 
For that we perform numerical analysis (again, both SOC 
and absorbing state experiments) using Eq.([6]). A di- 
rect integration of Eq.([5]) in one-dimension, using the re- 
cently introduced integration scheme for Langevin equa- 
tions with square-root noise [21| , produces the exponents 
reported in the first row of Table [V] (plots not shown) . 
All of them are compatible with those of the microscopic 
MZ model and the C-DP class (see Table H]). Changing 
the boundary conditions during the integration we im- 
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FIG. 3: Avalanche exponents for the anisotropic MZ sand- 
pile model in one dimension, averaged over 810 6 runs (system 
size L — 2 18 ). Up: spreading experiments (see Table 
Down: avalanche size (main plot) and time (inset) distribu- 
tions. Green lines mark DP scaling while red ones correspond 
to C-DP scaling. 
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0.28(5) 


0.21(5) 


1.14(5) 
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-0.39(5) 
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-0.01(2) 


0.50(2) 


1.39(5) 


1.64(5) 


1.98(2) 


0.51(2) 



TABLE V: Critical exponents for Eq. JB), with a reflective wall, 
Eq.lJBJl with an absorbing wall, and Eq.JB} with an anisotropic term. 



plement the reflecting or the absorbing wall. For the for- 
mer, we impose p{—x,i) — p{x,t), E(—x,t) — E(x,t), 
while for the absorbing walls p{x < 0,i) = 0, E(x < 
0, t) = 0. The measured exponents, performing avalanche 
and spreading experiments nearby a reflecting (absorb- 
ing) wall at (results not shown) are summarized in the 
second and third row of Table [Vj Again, the exponents 
coincide within error-bars with their corresponding C-DP 
counterparts and exclude DP scaling (see Table HT| . 

Finally, we have studied an anisotropic version of the 
equations by introducing a term proportional to V p(x, t) 
into both, the activity and the energy equations in 
Eq.®, obtaining again excellent agreement with the one- 
dimensional C-DP values (Table |LTT| . 

In summary, we have integrated numerically Eq.®, 



and implemented the necessary modifications (i.e. in- 
clude boundaries or anisotropy) to perform the tests de- 
scribed in the previous section. The obtained results 
are in excellent agreement with those for the microscopic 
model, confirming that (i) the Langevin equation derived 
in section II is representative of MZ model, and that (ii) 
the MZ model is in the C-DP class. 



VI. CONCLUSION AND DISCUSSION 

We have shed some light on the picture of universal- 
ity in stochastic sandpiles, by confirming that, indeed, 
they all share the same universal critical behavior. As 
hypothesized some years ago, their critical features are 
captured by the set of Langevin Eq.fl]), C-DP, describ- 
ing in a minimal way the phase transition into a mul- 
tiply degenerated absorbing state in the presence of a 
non-diffusive conserved field. 

We have shown that the Maslov-Zhang sandpile, be- 
lieved before to exhibit a different type of scaling (di- 
rected percolation like), is actually in the C-DP class, 
in agreement with the universality hypothesis. To reach 
this conclusion we have performed large scale simulations 
and introduced new numerical strategies to easily dis- 
criminate between DP and C-DP. In particular, we have 
benefited from the fact that the, otherwise very simi- 
lar, DP and C-DP classes behave in radically different 
ways both in the presence of walls and when anisotropy 
is switched on. 

We have also derived, in two different ways, an al- 
ternative set of Langevin equations, Eq.®, describing 
the Maslov-Zhang sandpile. This new set of equations 
is characterized by a different form of local energy diffu- 
sion (the corresponding current is proportional to energy 
gradients and not to activity gradients as is the case in 
Ecj-©)- By direct integration of the stochastic differen- 
tial set of Eq.©, we have shown that it describes the 
same universality class as Eq.((T]), i.e. C-DP, despite of 
the formal differences in their respective equations for the 
conserved field, hence leading to a coherent global picture 
for the universality of sandpiles. This result actually en- 
larges the C-DP universality class, allowing to embrace 
also different types of energy relaxation or redistribution 
dynamics. 

Our analyses have several general implications for the 
C-DP universality class: 

i) Reflecting walls arc not a relevant perturbation in 
this class: avalanche and spreading exponents measured 
in the vicinity of a reflecting wall coincide with their cor- 
responding bulk counterparts. The underlying reason for 
this remains to be well understood. 

ii) Absorbing walls are relevant ingredients and affect 
the corresponding surface critical behavior. In particular, 
avalanches and spreading experiments performed nearby 
the wall are characterized by exponents that differ from 
their bulk counterparts. 

iii) Anisotropy in space is also a relevant ingredient. 
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The corresponding critical behavior is described by the 
set of Langevin Eqs.([T]) (or, equivalently Eq.©) with an 
extra term Vp(x, t) in both equations. Contrarily to the 
isotropic case, the critical exponents in the anisotropic 
class are known exactly in any dimension. The results 
coincide with those of anisotropic interfaces in random 
media, confirming once again the equivalence between 
the absorbing-state and the interface pictures for SOC 
sandpiles [i~3| |. 

It would be highly desirable to have a working renor- 
malization group calculation allowing to put all the re- 



sults discussed here under a solid analytical ground. 
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Confirming and extending the hypothesis of universality in sandpiles 
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Stochastic sandpiles self-organize to an absorbing-state critical point with scaling behavior differ- 
ent from directed percolation (DP) and characterized by the presence of an additional conservation 
law. This is usually called C-DP or Manna universality class. There remains, however, an exception 
to this universality principle: a sandpile automaton introduced by Maslov and Zhang, which was 
claimed to be in the DP class despite of the existence of a conservation law. In this paper we show, 
by means of careful numerical simulations as well as by constructing and analyzing a field theory, 
that (contrarily to what previously thought) this sandpile is also in the C-DP or Manna class. This 
confirms the hypothesis of universality for stochastic sandpiles and gives rise to a fully coherent 
picture of self-organized criticality in systems with conservation. In passing, we obtain a number 
of results for the C-DP class and introduce a new strategy to easily discriminate between DP and 
C-DP scaling. 

PACS numbers: 05.50.+q,02.50.-r,64.60.Ht,05.70.Ln 



I. INTRODUCTION 

Aimed at shedding some light on the origin of scale 
invariance in many contexts in Nature, different mecha- 
nisms for self-organized criticality (SOC) were proposed 
during the last two decades following the seminal work by 
Per Bak and collaborators [l|, |2|, |3|]. Sandpiles, ricepiles, 
and earthquake toy-models become paradigmatic exam- 
ples capturing the essence of self-organization to scale- 
invariant (critical) behavior without apparently requiring 
the fine tuning of parameters 0, d, 0] . 

Sandpiles played a central role in the development of 
this field P, H, Q • They are metaphors of real systems 
(as earthquakes, snow avalanches, stick-slip phenomena, 
etc.) in which some type of stress or energy is accumu- 
lated at some slow timescale and relaxed in a much faster 
way. In sandpiles, grains are slowly added until eventu- 
ally they relax if a local instability threshold is overcome; 
then, they are transmitted to neighboring sites which, on 
their turn, may become unstable and relax, generating 
avalanches of activity. Considering open boundaries to 
allow for energy balance, a steady state with power-law 
distributed avalanches is eventually reached. 

In order to rationalize sandpiles in particular and SOC 
in general and to understand their critical properties, it 
was proposed @, [1, Q to look at them as systems with 
many absorbing states [l(J [H| • The underlying idea is 
that, in the absence of external driving, sandpile mod- 
els get eventually trapped into stable configurations from 
which they cannot escape, i.e. absorbing states [To] . The 
concept of fixed energy sandpiles (FES) was introduced 
to make this connection more explicit. FES sandpiles 
share the microscopic rules with their standard (slowly 
driven and dissipative) counterparts, but with no driving 
nor dissipation. In this way, the total amount of sand 
or energy becomes a conserved quantity acting as a con- 
trol parameter. Calling E c the average density of energy 



of a standard sandpile in its stationary (self-organized) 
critical state, it has been shown that the corresponding 
fixed-energy sandpile exhibits a transition from an active 
phase to an absorbing one at, precisely, E c , while it is ab- 
sorbing (active) below (above) E c . It can be argued that 
slow driving and dissipation acting together at infinitely 
separated time-scales constitute a mechanism able to pin 
a generic system with absorbing states and a conservation 
law to its critical point [H, @, 0, H, S, [HJ • 

Usingthe relation with transitions into absorbing 
states |l3[ , a Langevin equation describing FES stochas- 
tic sandpiles was proposed [1, Q in the spirit of Ho- 
henberg and Halperin [l4j. The Langevin equation 
is similar to the well known directed-percolation (DP) 
Langevin equation describing generic systems with ab- 
sorbing states [HI EH but it is coupled linearly to a con- 
served non-diffusive energy field, representing the conser- 
vation of energy grains in the sandpile dynamics d, [§] : 

d t p(x,t) = ap - bp 2 + ojpE(x,t) + X7 2 p + a^/p?](x,t) 
d t E(x,t) = DV 2 p, (1) 

where D, a, b and lo are constants, p(x,t) and E(x,t) 
are the activity and the energy field, respectively and 
r\ is a zero-mean Gaussian white noise. The universality 
class described by this Langevin equation, which includes 
sandpiles as well as all other systems with many absorb- 
ing states and a non-diffusive conserved field, is usually 
called Manna or C-DP class 

US, HE HI- As 

an impor- 
tant side note, let us remark that the original Bak-Tang- 
Wiesenfeld model, being deterministic, has many other 
conservation laws (toppling invariants) and is therefore 
not described by the present stochastic theory 0, EH . 

Incidentally, even if it has been clearly established that 
DP and C-DP constitute two different universality classes 
[l9L I20L l2lj , most of their universal features (critical ex- 
ponents, moment ratios, scaling functions,...) are very 
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similar, making it difficult to discriminate numerically 
between both classes in any spatial dimension. 

Despite of the fact that the critical behavior of stochas- 
tic sandpiles is accepted to be universal and described by 
Eq.([T]), there are a few sandpile models which, despite of 
being conservative, have been argued to exhibit a differ- 
ent type of critical behavior [22, 2j| , namely DP, violating 
the conjecture of universality: 

• the Mohanty-Dhar (MD) sandpile 22], in which 
grains have some probability to remain stable even 
if they are above the instability threshold. It was 
argued, based on numerical simulations and on ex- 
act results for an anisotropic version of it, that the 
(isotropic) MD sandpile should exhibit DP behav- 
ior, 

• the Maslov-Zhang (MZ) sandpile (23|, in which not 
only the grains of unstable sites are redistributed 
among nearest neighbors. Instead, all the grains 
in the local neighborhood of unstable sites are ran- 
domly redistributed or "reshuffled" . Various ver- 
sions of the model (called "charitable" , "neutral" 
and "greedy") were defined depending on the bias 
in the local redistribution rule ( "charitable" if the 
central site receives less than each n.n., "neutral" if 
all the sites in the neighborhood are treated alike, 
and "greedy" otherwise). On the basis of Monte 
Carlo simulations this sandpile (in its neutral ver- 
sion) was argued to be DP-like. 

Aimed at clarifying this puzzling situation, in a recent 
work we provided strong numerical and analytical ev- 
idence that, contrarily to what previously thought, the 
(isotropic) MD sandpile is actually in the C-DP class [24| 
(see [26[ for a discrepant viewpoint). To further substan- 
tiate our claim, in [25| we proposed a strategy to easily 
discriminate between DP and C-DP consisting in intro- 
ducing a wall (either absorbing or reflecting); systems in 
the DP and in the C-DP classes behave in (qualitative 
and quantitatively) very different ways in the presence of 
walls, providing an easy criterion to discriminate between 
both classes. 

Our goal in the present paper is to scrutinize the last 
remaining piece in the puzzle, i.e. the MZ model, by em- 
ploying extensive numerical simulations as well as field 
theoretical considerations. For the numerics, we take ad- 
vantage of the previously introduced discrimination crite- 
rion and, also, present another method to easily discrim- 
inate between DP and C-DP, based in the introduction 
of anisotropy in one spatial direction. 

The analyses presented on what follows, show in a 
clean-cut way that the MZ model, in any of the three ver- 
sions above, is actually in the C-DP class. This leaves no 
dangling end in the sandpile universality picture and, as a 
byproduct, confirms the general validity of the absorbing- 
state approach to rationalize SOC. New results for the 
C-DP class are also obtained. 



II. THE MASLOV-ZHANG SANDPILE 

The (neutral) Maslov-Zhang (MZ) sandpile is de- 
fined as follows: 

1. Driving: an input energy, 5E < 1, is added to the 
central (or to a randomly chosen) site, i, of a d- 
dimensional lattice, and the site is declared active. 

2. Relaxation: energy is locally redistributed ( "reshuf- 
fled") between the active site and its nearest- 
neighbors, according to: 

2d+l 

E i = E ( 2 ) 

where are uniformly distributed (rj € [0, 1]) ran- 
dom variables, and the sums are performed over 
the site i and its 2c? nearest neighbors. This rule 
needs to be slightly modified for the "charitable" 
and "greedy" versions of the model that we will 
not explore in detail here. 

3. Activation: each of the sites involved in the reshuf- 
fling is declared active with a probability given 
by its own energy, triggering the generation of 
avalanches of activity. 

4. Avalanches: new active sites are added to a list and 
relaxed in a sequential way. 

5. Dissipation: energy arriving at the open borders is 
removed from the system. 

6. Avalanches proceed until all activity ceases, and 
then a new external input is added. Eventually a 
critical stationary state is reached. 

Monte Carlo simulations by Maslov and Zhang re- 
vealed exponents compatible with those of directed per- 
colation in two dimensions and above, while in one di- 
mension some anomalies in the scaling were reported 
[23l ]. Later, simulations of the FES counterpart of the 
MZ sandpile led to very similar results [8|. 

Before entering a more careful numerical analysis of the 
MZ sandpile, we take a detour to construct explicitely a 
Langevin equation for such a model. This will allow us 
to better understand what the main relevant ingredients 
of the MZ model are and in what sense it differs from 
other C-DP models and from Eq.{T]). 

III. A LANGEVIN EQUATION FOR THE 
MASLOV-ZHANG SANDPILE 

The main difference between the Maslov-Zhang cellu- 
lar automaton and other sandpiles in the C-DP univer- 
sality class is that, while in the C-DP class as the Manna 
model, the only energy redistributed by the dynamics (by 
topplings) is that accumulated in active sites, in the MZ 
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dynamics both the energy of active sites as well as that 
of its nearest neighbors is redistributed. This leads to a 
more severe local redistribution of energy, that we quan- 
tify in what follows in terms of a new set of Langevin 
equations, first in a phcnomcnological way and then by 
deriving it from the microscopic rules. 



A. Phenomenological Langevin equation 

Diffusion of energy occurs in the C-DP class by means 
of activity relaxation. At a mesoscopic level, this implies 
that the rate of change of the energy density E(x, t) at a 
given position is proportional to minus the divergence of 
a current, dtE(x,t) = — V.j(x,t), where j(x,t) is given 
by the gradient of the activity field, 



j(x,t) = -DVp(x,t) 



(3) 



leading to d t E(x, t) = D\7 2 p(x, t) in Eq.©. 

Instead, in the MZ model, changes of energy are con- 
trolled by the reshuffling rule, i.e. energy does not need 
to be at an active site (but in its neighborhood) to be 
redistributed. Hence, at a mesoscopic level the current 
is given by gradients of the energy in the presence of 
non-vanishing activity. More specifically, dtE(x, t) = 
— V.j(cc,i) where now 



j(x,t)=-D VE(x,t), 



(4) 



and the diffusion D is not a constant but a functional pro- 
portional to p(x,t), D(p(x,t)) = Dp(x,t) capturing the 
requirement that local reshuffling of energy only occurs 
in the presence of activity. This enforces the absorbing 
state condition that dynamics ceases if p = 0. Finally, 



d t E{x,t) = D V.[p{x,t)VE(x,t)]. 



(5) 



On the other hand, the equation for the activity is not 
expected to change in any relevant way, so one obtains 

d t p(x,t) = ap - bp 2 + V 2 p + ujpE(x,t) + a^/pn(x,t) 
d t E(x,t) = DV.[p(x,t)VE{x,t)} (6) 

representing the MZ dynamics at a mesoscopic level. 
This is to be compared with Eq. (JTJ) . 



B. Microscopic derivation of the Langevin equation 

To gain more confidence on the phenomenological set 
of equations, Eq. ([6]) we present here an explicit micro- 
scopic derivation. We consider a parallel version of the 
MZ model in which all active sites are relaxed at every 
time step. To do so, we assume that each site uses a 
fraction l/(2d+ 1) of its total energy for eventual redis- 
tributions with each of the sites in its local neighborhood. 



At mean-field level, the energy evolves according to 

2d+l 

- — K, , 

2d 



E, 



i,t + l 



1 2d+l 

—E itt £ [i - e( Pj , t ) 



2d+l 



2d+l 



2d 



E (^) E 



(7) 



i=i 



k=l 



where j (respectively k) runs over the nearest neighbors 
of i (j) and <d is the Heaviside step function (6(z) = 
if z < 0). For any inactive site in the local neighborhood 
(i.e. pj = 0), the central site does not redistribute the 
corresponding fraction of its energy (first term in Eq.©). 
For any active site in the local neighborhood (i.e. pj > 
0) the central site receives on average a corresponding 
fraction of energy (second term). Eq.© is a mean- field 
equation to which fluctuations should be added to have 
a more detailed description. Such fluctuations appear 
as a (conserved) noise for the resulting energy equation 
and can be easily argued to constitute a higher order, 
irrelevant, correction. 

Regularizing the Heaviside step function in Eq.© by 
means of a hyperbolic-tangent and expanding it in power 
series up to first order in p [24|j we obtain 
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Introducing the d-dimcnsional discrete Laplacian V 2 Ei = 
y^j—i (Ej — Ei) the first term can be rewritten as: 



E, 



2d+l 



2d 



Ei 



E, 



2d+l 



2d 



I E Pj 



(9) 



2d+ 1 

Similarly, the second term can be expressed as 



(2d 



2d+l /2d+l \ 

^E^(E^J=^ + ^TT)^(^) 

(10) 

Putting these two contributions together and reorganiz- 
ing the discrete derivatives, one obtains: 



Ei,t+i = Eij + 



1 



2d+ 1 



pi,tV 2 E i>t 



2d+l 



Vpi,t S7Ei. 



Eit 



{2d+ jV 2 (P^E ht ) = 



2d+l 



(ii) 
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which in the continuous time limit becomes 
2 _ , _„ , 1 



d t Ei = 



2d+l 



{2d+l)< 



V 2 U t V 2 E ht ) 



(12) 

The second term in this last equation can be argued to be 
naively irrelevant in the renormalization group sense (as 
it includes higher order derivatives) and hence dropped 
out. Identifying D with 2/{2d + 1) we recover the phe- 
nomenological Eq.©, as the leading contribution. 

As said above, this is to be compared with the stan- 
dard Langevin Eq.© for the C-DP class, in which the 
flux of energy is given by the gradient of the activity 
itself. The following questions pop up naturally: Does 
Eq. © lead to a critical behavior different from that of 
Eq.©? Could this be the reason why the MZ sandpilc 
was claimed to differ from the C-DP class? Is this new 
form of conserved energy dynamics irrelevant at the DP 
fixed point, supporting the MZ sandpile to be DP like? 

To properly answer these questions one should resort 
to a full renormalization group calculation. Given that 
even for the C-DP class this has proven to be a, still 
not satisfactorily accomplished, difficult task [13] , we will 
leave aside such a strategy here. Instead, in section IV 
we will give an answer to these question by means of 
computational studies of the microscopic MZ model as 
well as of its equivalent Langevin Eq.©. 



C. Naive power counting 

A power counting analysis is not helpful in elucidating 
the relevancy of the energy-diffusion term, Eq. © . Actu- 
ally, as there is a lineal dependence on the energy field in 
both the r.h.s. and the l.h.s. of Eq.©, there is no way 
to extract the energy field naive dimension nor, hence, to 
make any statement about the relevancy of the coupling 
term wp(x,t)E(x,t) at the DP fixed point. 

On the other hand, it is easy to cast Eq.© into a 
generating functional following standard procedures, to 
set the basis of a perturbative expansion. However, one 
soon realizes that technical difficulties similar to those 
encountered for the analysis of Eq. ([6]) (including the pres- 
ence of generically singular propagators H,[13]) show up, 
hindering the perturbative calculation. It is our believe 
that some type of non-perturbative technique, or non- 
conventional perturbative expansion, is required to elu- 
cidate the renormalization group fixed point of this type 
of problems. 

The analytical understanding, at a field theory level of 
C-DP as well as the MZ-Langevin equation remain open 
challenging tasks. 



its asymptotic (long time and large system size) prop- 
erties. We report on two different types of numerical 
experiments. 

• "SOC" or avalanche experiments: 

By iterating slow addition of grains in the sandpilc 
with open boundaries, the system self-organizes 
to a state with average energy E c . Then, the 
avalanche size distribution, P(s) and the avalanche 
time distribution P t (t) , can be estimated and their 
corresponding exponents, r and Tt [28| . measured. 

• Absorbing state experiments: 

At the stationary state, we perform spreading ex- 
periments from a localized seed, and we measure: 
i) the mean quadratic distance to the initial seed 
R 2 t Zapr in active runs, ii) the average number of 
active sites as a function of time, N(i) ~ t v and iii) 
the surviving probability up to time t, P s ~ t~ s . 
We also study the decay at criticality of a homo- 
geneous initial activity (p(t) ~ t~ e ) in the fixed 
energy case [l(| HH • 

Scaling laws relating avalanche to spreading exponents 
where described systematically in [2jj; two of them are: 
r = (l + r/ + 2<5)/(l + r/ + <5), T t = 1 + S. We measure the 
exponents independently and use scaling laws as a check 
for consistency. 

We simulate the MZ automaton in one-dimensional lat- 
tices up to size L = 2 15 . The stationary critical energy 
density is E c = 0.4928(2) and, contrarily to what re- 
ported in [23|, we do observe clean scaling at criticality, 
even if it emerges only after significatively long transients 
(results not shown) justifying why smaller scale simula- 
tions can lead to erroneous conclusions. 

The resulting critical exponents are gathered together 
in Table Q] They are closer in all cases to the C-DP 
values than to DP ones. The exponents measured ex- 
plicitely in [H are a (J dtN(t)/P s (t) ~ t a ) and the 
fractal dimension Df, which take also very similar val- 
ues in both cases (a = 1.47(1), D f = 2.32(1) for DP, 
a = 1.52(1), D f = 2.36(1) for C-DP). Still, given the 
similarity between the numerical values in both classes, 
it is not save to extract a definitive conclusion from these 
results. Higher numerical precision would be required to 
produce fully convincing evidence. As said above, the 
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DP 


0.313(1) 


0.159(1) 


1.108(1) 


1.159(1) 


1.265(1) 


0.159(1) 


C-DP 


0.350(5) 


0.170(5) 


1.11(2) 


1.17(2) 


1.39(1) 


0.125(1) 


MZ 


0.32(5) 


0.20(5) 


1.13(5) 


1.20(5) 


1.40(5) 


0.13(1) 



IV. MONTE CARLO SIMULATIONS OF THE 
MASLOV-ZHANG SANDPILE 

We have performed extensive Monte Carlo simulations 
of the MZ sandpile (and variations of it) and scrutinized 



TABLE I: Critical exponents for DP, C-DP, and the MZ sandpile 
in one dimension. DP and C-DP values taken from 20, 25, 281. 



numerical values of DP and C-DP exponents are closer 
and closer as the dimensionality is increased (they co- 
incide above d c = 4). Therefore, performing numerical 
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simulations, as the ones above, to discriminate between 
both classes by using larger and larger times and system 
sizes in d > 2 , is not a clever idea. Instead, it is advis- 
able to use/devise more effective numerical strategies to 
discriminate between DP and C-DP in a simple, efficient, 
and numerically inexpensive way. For this we use: 

• the method devised in [25| consisting in analyzing 
how the system responds to the presence of a wall, 
or 
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- Ps(t) 




DP 




- C-DP 





• a new strategy, which exploits the fact that systems 
in these two classes react in remarkably different 
ways to the introduction of anisotropy in space. 

Both of these strategies allow to obtain clean-cut re- 
sults, as shown in the forthcoming two subsections. 

A. Boundary Driven Experiments 

The influence of walls in systems in the DP class has 
been profusely analyzed in the literature |29j ] . In par- 
ticular, it is well known that, if spreading (and SOC) 
experiments are performed nearby a wall, the surviving 
probability is significatively affected, and avalanche and 
spreading exponents change in a non-trivial way with re- 
spect to their bulk counterparts [29J. It is also known 
that in the DP class, both reflecting and absorbing walls 
lead to a common type of universal "surface critical be- 
havior", that we call surface directed percolation (SDP), 
characterized (in one-dimension) by the exponents shown 
in Table II. 

In contrast, the effect of walls in C-DP systems has 
been studied only recently [25]. Contrarily to the DP 
case, absorbing and reflecting walls induce different types 
of surface critical behavior. As illustrated in Table II, all 
spreading and avalanche exponents take distinct values 
for an absorbing and for a reflecting wall. Furthermore, 
the numerical differences between the exponents for ei- 
ther type of wall with respect to their corresponding SDP 
counterparts are very large, allowing for easy numerical 
discrimination [25]. Finally, in the C-DP class, the ex- 
ponents in the presence of a reflecting wall coincide with 
their bulk counterparts 25]. These features imply that, 
by introducing a wall in a given system with absorbing 
states, it becomes straightforward to distinguish if it is 
in the DP or in the C-DP class, with moderate compu- 
tational cost. 

Following this strategy, we simulated the one- 
dimensional MZ sandpile, as defined above, in the pres- 
ence of both reflecting and absorbing walls. In both cases 
a wall is introduced at the origin (position i = 0), and the 
sandpile is studied in the positive half lattice. In the re- 
flecting case, the energy that should go after reshuffling, 
to the leftmost site, at i = 0, (whose energy is fixed to 
zero) is instead added to its closest nearest neighbor to 
the right, i = 1. On the other hand, the absorbing condi- 
tion is imposed by fixing the energy of the leftmost site 
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FIG. 1: (Color online) Avalanche exponents for the one- 
dimensional MZ sandpile in the presence of a reflecting wall. 
Up: spreading experiments; N(t) is the total number of ac- 
tive sites and Ps(t) is the surviving probability; system size 
L — 2 15 number of runs 500 . Down: avalanche size (main 
plot) and time (inset) distributions; system size L — 2 12 num- 
ber of runs 10 5 . Green (red) lines mark DP (C-DP) scaling 
(see numerical values in Table [TTJ . 



to zero after every iteration of the microscopic sandpile 
rules, i.e. by removing from the system at every iteration 
all the energy received by the leftmost site. 
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0.045(2) 


0.425(2) 
0.426(2) 
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1.426(2) 
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C-DP re/ 
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0.35(1) 
-0.33(2) 


0.16(1) 
0.85(2) 


1.11(2) 
1.56(2) 


1.15(2) 
1.81(2) 


1.41(1) 
1.43(2) 


MZ re/ 
MZ abs 


0.37(5) 
-0.36(5) 


0.15(5) 
0.84(5) 


1.09(5) 
1.57(5) 


1.14(5) 
1.84(5) 


1.30(5) 
1.25(5) 



TABLE II: One-dimensional critical exponents for DP and C-DP 
without walls 20, 28] and in the presence of absorbing and reflect- 
ing walls [25ll , Values in rows 1 (DP re ^) and 2 (DP ai)s ) coincide 
within error-bars. Note also, that values in row 3 ( C-DP re ^ ) coin- 
cide with those for C-DP (table [IJ. Results for the MZ sandpile in 
the presence of reflecting/absorbing walls are reported in the last 
two rows. 

Fig.{T]) and Fig. ([2]) show the results of simulations per- 
formed in lattices of system size L = 2 15 , averaging over 
up to 8 • 10 7 runs. The corresponding exponents are sum- 
marized in Table II. All of them coincide within numerical 
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FIG. 2: (Color online) As in Fig. [T] but for an absorbing 
wall (numerical values summarized in Table [II}. System size 
L = 2 15 and 8 ■ 10 7 runs. 

accuracy with the expected values for the C-DP class in 
the presence of reflecting or absorbing walls respectively, 
and differ significatively from those of SDP. For example, 
in the presence of a reflecting (respectively, absorbing) 
wall, the measured value of r\ is 0.37 (—0.37) in good 
agreement with the C-DP expectation, r\ — 0.35 (—0.33) 
and in blatant disagreement with the corresponding DP 
value, r\ « 0.046 (« 0.045), which is one order of mag- 
nitude smaller (and of opposite sign in the case of an 
absorbing wall). Similar large differences are measured 
for all the exponents (see table II). Note also that, as is 
the case in the C-DP class [25| , the exponents in the pres- 
ence of a reflecting wall coincide within error-bars with 
their bulk counterparts. In conclusion, studying the in- 
fluence of walls we conclude that the one-dimensional MZ 
sandpile exhibits C-DP scaling. 



B. Anisotropic Experiments 

It is well known that systems in the DP class are invari- 
ant under Galilean transformations: if particles have a 
tendency to move anisotropically in one preferred spatial 
direction, that does not alter the critical properties [l(J. 
The presence of any degree of anisotropy in DP-like sys- 
tems is an irrelevant trait, or in other words, anisotropic 
DP (A- DP) is just DP. 

The role of anisotropy in sandpiles has also been pro- 
fusely studied after the pioneering exact solution by Dhar 



and Ramaswamy [3(| of the totally anisotropic or "di- 
rected" counterpart of the Bak-Tang-Wiescnfcld sand- 
pile. Anisotropic stochastic sandpiles have also been 
studied using general principles [3l| and through in- 
terfacial representations (32J. The conclusion is that 
all anisotropic sandpiles, as long as they are stochastic 
[33| . belong to the same universality class, that we call 
anisotropic C-DP (A-C-DP) [34(. The critical exponents 
of models in this class, where first measured numerically 
(33 | . and then exactly calculated in any dimension [35J ] 
(see table HTT1 and table ITVT). 
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TABLE III: One dimensional critical exponents for D P |20t |28(| . 
C-DP with a preferred direction (analytical results from [351]). and 
the anisotropic MZ model. 
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TABLE IV: Two-dimensional critical exponents for anisotropic 
models: A-DP (i.e. DP), A-C-DP (analytical results), and 
anisotropic-MZ. 



The strategy to be used is straightforward: take the 
MZ sandpile model and switch on anisotropy; if the 
isotropic model was in the DP class, anisotropy should be 
an irrelevant ingredient and the anisotropic counterpart 
should also be DP like. If, instead, the isotropic model is 
in the C-DP class, then anisotropy is a relevant ingredi- 
ent and critical exponent change from C-DP to A-C-DP 
values. 

The simplest way to define an anisotropic MZ (A-MZ) 
model is by fixing one of the rj in Eq. , say the one to 
the right, to its maximum possible value, r ■ = 1, and let- 
ting the others rj to take randomly distributed values in 
[0, 1]. This generates an overall energy flow towards the 
preferred direction (to the right, in this case). Anisotropy 
can be introduced in other ways, including full-anisotropy 
or directness, but this does not affect our conclusions in 
any significant way. 

Fig. [3] and table Hill show our main results for the one- 
dimensional MZ model with anisotropy. Both avalanche 
and spreading exponents are very different from their 
isotropic counterparts. They also differ notoriously from 
DP values, but coincide within error-bars with the ex- 
pected values for the A-C-DP class. The same conclusion 
holds in two dimensions (see Table IIV[) . In this way, as 
the anisotropic MZ model belongs to A-C-DP class the 
original, isotropic, MZ sandpile model can be safely con- 
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FIG. 3: (Color online) Avalanche exponents for the 
anisotropic MZ sandpile model in one dimension, averaged 
over 8 • 10 6 runs (system size L — 2 18 ). Up: spreading exper- 
iments (see Table ITTT]) . Down: avalanche size (main plot) and 
time (inset) distributions. Green lines mark DP scaling while 
red ones correspond to C-DP scaling. 



eluded to be in the C-DP universality class, confirming 
the result above. 



V. NUMERICAL INTEGRATION OF THE MZ 
LANGEVIN EQUATION 

In this section we verify that Langevin Eq.© is a 
sound description of the MZ model and that, despite of 
its different form, it behaves asymptotically as Eq.©. 
For that we perform numerical analysis (again, both SOC 
and absorbing state experiments) using Eq.©. A di- 
rect integration of Eq.© in one-dimension, using the re- 
cently introduced integration scheme for Langevin equa- 
tions with square-root noise [2l| , produces the exponents 
reported in the first row of Table |V] (plots not shown) . 
All of them are compatible with those of the microscopic 
MZ model and the C-DP class (see Table IJ). Changing 
the boundary conditions during the integration we im- 
plement the reflecting or the absorbing wall. For the for- 
mer, we impose p{—x,i) — p{x,t), E(—x,t) — E(x,t), 
while for the absorbing walls p(x < 0,t) — 0, E(x < 
0, t) = 0. The measured exponents, performing avalanche 
and spreading experiments nearby a reflecting (absorb- 
ing) wall at (results not shown) are summarized in the 
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TABLE V: Critical exponents for Eq.JBJ, with a reflective wall, 
Eq. (O with an absorbing wall, and Eq. J6} with an anisotropic term. 



second and third row of Table [V] Again, the exponents 
coincide within error-bars with their corresponding C-DP 
counterparts and exclude DP scaling (see Table [TTJ) . 

Finally, we have studied an anisotropic version of the 
equations by introducing a term proportional to V p(x,t) 
into both, the activity and the energy equations in 
Eq.®, obtaining again excellent agreement with the one- 
dimensional C-DP values (Table [ill]) . 

In summary, we have integrated numerically Eq.©, 
and implemented the necessary modifications (i.e. in- 
clude boundaries or anisotropy) to perform the tests de- 
scribed in the previous section. The obtained results 
are in excellent agreement with those for the microscopic 
model, confirming that (i) the Langevin equation derived 
in section II is representative of MZ model, and that (ii) 
the MZ model is in the C-DP class. 



VI. CONCLUSION AND DISCUSSION 

We have shed some light on the picture of universal- 
ity in stochastic sandpiles, by confirming that, indeed, 
they all share the same universal critical behavior. As 
hypothesized some years ago, their critical features are 
captured by the set of Langevin Eq.JT]), C-DP, describ- 
ing in a minimal way the phase transition into a mul- 
tiply degenerated absorbing state in the presence of a 
non-diffusive conserved field. 

We have shown that the Maslov-Zhang sandpile, be- 
lieved before to exhibit a different type of scaling (di- 
rected percolation like), is actually in the C-DP class, 
in agreement with the universality hypothesis. To reach 
this conclusion we have performed large scale simulations 
and introduced new numerical strategies to easily dis- 
criminate between DP and C-DP. In particular, we have 
benefited from the fact that the, otherwise very simi- 
lar, DP and C-DP classes behave in radically different 
ways both in the presence of walls and when anisotropy 
is switched on. 

We have also derived, in two different ways, an al- 
ternative set of Langevin equations, Eq.®, describing 
the Maslov-Zhang sandpile. This new set of equations 
is characterized by a different form of local energy diffu- 
sion (the corresponding current is proportional to energy 
gradients and not to activity gradients as is the case in 
Eq.([T])). By direct integration of the stochastic differen- 
tial set of Eq.©, we have shown that it describes the 
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same universality class as Eq.([T]), i.e. C-DP, despite of 
the formal differences in their respective equations for the 
conserved field, hence leading to a coherent global picture 
for the universality of sandpiles. This result actually en- 
larges the C-DP universality class, allowing to embrace 
also different types of energy relaxation or redistribution 
dynamics, which includes also "charitable" versions of 
the MZ model. Instead, the "greedy" version, charac- 
terized by "anti-diffusion" (i.e. energy accumulates on 
active sites) is expected to be highly anomalous. 

Our analyses have several general implications for the 
C-DP universality class: 

i) Reflecting walls are not a relevant perturbation in 
this class: avalanche and spreading exponents measured 
in the vicinity of a reflecting wall coincide with their cor- 
responding bulk counterparts. The underlying reason for 
this remains to be well understood. 

ii) Absorbing walls are relevant ingredients and affect 
the corresponding surface critical behavior. In particular, 
avalanches and spreading experiments performed nearby 
the wall are characterized by exponents that differ from 
their bulk counterparts. 

iii) Anisotropy in space is also a relevant ingredient. 



The corresponding critical behavior is described by the 
set of Langevin Eqs.([T]) (or, equivalently Eq.©) with an 
extra term V p{x,t) in both equations. Contrarily to the 
isotropic case, the critical exponents in the anisotropic 
class are known exactly in any dimension. The results 
coincide with those of anisotropic interfaces in random 
media, confirming once again the equivalence between 
the absorbing-state and the interface pictures for SOC 
sandpiles [Ta j". 

It would be highly desirable to have a working renor- 
malization group calculation allowing to put all the re- 
sults discussed here under a solid analytical ground. 
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